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Abstract. We use a three-dimensional lattice Boltzmann model to in- 
vestigate the spreading of mesoscale droplets on homogeneous and het- 
erogeneous surfaces. On a homogeneous substrate the base radius of the 
droplet grows with time as t ' 28 for a range of viscosities and surface 
tensions. The time evolutions collapse onto a single curve as a function 
of a dimensionless time. On a surface comprising of alternate hydropho- 
bic and hydrophilic stripes the wetting velocity is anisotropic and the 
equilibrium shape of the droplet reflects the wetting properties of the 
underlying substrate. 



1 Introduction 

Wetting processes, such as the spreading of a droplet over a surface, have at- 
tracted the attention of scientists for a long time [1] . A great deal is understood 
about the wetting behaviour of equilibrium droplets. However less is known 
about the dynamics of these systems, a problem of considerable industrial rel- 
evance with the advent of ink-jet printing. The droplets involved in printing 
typically have length scales of microns. Experimental work on such mesoscopic 
droplets is difficult and expensive because of the length- and time-scales involved. 
Therefore there is a need for numerical modelling both to investigate the physics 
and to help design and interpret the experiments. 

Lattice Boltzmann models are a class of numerical techniques ideally suited 
to probing the behaviour of fluids on mesoscopic length scales [2]. Several lat- 
tice Boltzmann algorithms for a liquid-gas system have been reported in the 
literature [3-5]. They solve the Navier-Stokes equations of fluid flow but also 
input thermodynamic information, typically either as a free energy or as ef- 
fective microscopic interactions. They have proved successful in modelling such 
diverse problems as fluid flows in complex geometries [6], two-phase models [3, 
4] , hydrodynamic phase ordering [7] and sediment transport in a fluid [8] . 

Here we show that it is possible to use a lattice Boltzmann approach to model 
the spreading of mesoscale droplets and, in particular, to illustrate how a droplet 
spreads on a substrate comprising of hydrophilic and hydrophobic stripes. 
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We consider a one-component, two-phase fluid and use the free energy model 
originally described by Swift et al. [3] with a correction to ensure Galilean invari- 
ance [9] . The advantage of this approach for the wetting problem is that it allows 
us to tune equilibrium thermodynamic properties such as the surface tension or 
static contact angle to agree with analytic predictions. Thus it is rather easy to 
control the wetting properties of the substrate. Three dimensional simulations of 
spreading on smooth and rough substrates have previously been reported in [10] 
using a different lattice Boltzmann algorithm. 

The paper is organised as follows. First we summarise the main features of the 
lattice Boltzmann approach. The model is validated by showing the consistency 
of the measured equilibrium contact angle with that predicted by Young's law 
and by measuring the base radius of the spreading droplet as a function of time 
obtaining, as expected, a power law growth. We show that when the reduced 
base radius is plotted as a function of reduced time the data fall on a universal 
curve for several values of surface tension and viscosity. 

We then consider spreading on a heterogeneous substrate consisting of alter- 
nate hydrophobic and hydrophilic stripes. We find that the spreading velocity 
is anisotropic and that the final droplet shape reflects the wetting properties of 
the underlying substrate. Finally, a conclusion suggests extensions to the work 
presented here. 

2 Simulating spreading 

2.1 The lattice Boltzmann model 

The lattice Boltzmann approach solves the Navier-Stokcs equations by following 
the evolution of partial distribution functions /, on a regular, d-dimensional 
lattice formed of sites r. The label i denotes velocity directions and runs between 
and z. DdQz+1 is a standard lattice topology classification. The D3Q15 lattice 
topology we use here has the following velocity vectors v,: (0, 0, 0), (±1, ±1, ±1), 
(±1,0,0), (0,±1,0), (0,0, ±1) in lattice units. 
The lattice Boltzmann dynamics are given by 



where At is the time step of the simulation, r the relaxation time and f^ q the 
equilibrium distribution function which is a function of the density n — J2i=o fi 
and the fluid velocity u defined through the relation rau = J2t=o /* Vi - 
The relaxation time tunes the kinematic viscosity as 
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where Ar is the lattice spacing and C2 and C4 are coefficients related to the 
topology of the lattice. These are equal to 3 and 1 respectively when one considers 
a D3Q15 lattice (see [11] for more details). 



It can be shown [3] that equation (1) reproduces the Navier-Stokes equations 
of a non-ideal gas if the local equilibrium functions are chosen as 

f- q = A a + B a u a v, ia + G CT u 2 + D a u a Uf3V ia Vif3 + G aaf 3V ia Vip, i > 0, 

/o e, =»-E/r (3) 
i=i 

where Einstein notation is understood for the Cartesian labels a and j3 (i.e. 
Vi a u a = ^2 a Vi a u a ) and where a labels velocities of different magnitude. 

The coefficients A a , B ai C a , D a and G CT are chosen so as to satisfy the 
relations 
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where P a p is the pressure tensor, c is defined to be Ay/ At and the last term of 
the third expression in equation (4) is included to ensure Galilean invariance. 

Considering a D3Q15 lattice and a square-gradient approximation to the 
interface free energy {n{d a n) 2 /2) [3], a possible choice of the coefficients is [12] 

A a — -^r (pb - ^(d a n) 2 - nnd aa n + vu a d a nj , 
_ w a n w a n _ 3w a n 
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Gi 77 = — ^ (n(d 7 n) 2 + 2vu 1 d 1 n) , G2 77 = 0, 

G2 7 a = (K(9 7 n)(9an) + v(u 7 dsn + usd 7 n)) (5) 

where toi = 1/3, u> 2 = 1/24, k is a parameter related to the surface tension and 
Pb is the pressure in the bulk which is defined below. One can show [13] that the 
pressure tensor can be written as 

P a p = (pb - |(3 7 n) 2 - Knd 71 nJ 5 a p + K(d a n){dpn). (6) 



2.2 Wetting boundary conditions 



In this paper, we will focus our attention on flat substrates normal to the z 
direction. The derivatives in that direction should then be handled in such a 



way that the wetting properties of the substrate can be controlled. A boundary 
condition can be established using the Cahn model [14]. He proposed adding 
an additional surface free energy <?c(n s ) = fa — <j)in s + ■ ■ ■ at the solid surface 
where n s is the density at the surface. Neglecting the second order terms of ^ c (n) 
and minimizing <Pf, + ^> c (where is the free energy in the bulk), a boundary 
condition valid at z = emerges [15] 

d z n = -^. (7) 

K 

Equation (7) is imposed on the substrate sites to implement the Cahn model in 
the lattice Boltzmann approach. Details are given in [15]. 

The Cahn model can be used to relate fa to 9 the contact angle defined as 
the angle between the tangent plane to the droplet and the substrate. Within 
the Cahn model the surface tension at the interfaces is given by [1] 
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where W(n, T) is the excess free energy, <r; s , a sg , <j s i are the surface tensions 
at the liquid-gas, solid-gas and solid- liquid interface respectively and n s , ni, 
n g are the densities at the substrate, of the liquid phase and of the gas phase 
respectively. Young's law [16] gives a relation between the static contact angle 
and the surface tensions of the three phases 

aig cos 6 — a sg - a sl . (9) 

A convenient choice of bulk pressure is [15] 

Pb = Pc{v P + l) 2 (3fp - 2v p + 1 - 2/?t p ) (10) 

where v v = (n — n c )/n c , t p = (T c — T)/T c and p c = 1/8, n c = 3.5 and T c = 4/7 
are the critical pressure, density and temperature respectively and (3 is a constant 
typically equal to 0.1. The excess free energy then becomes [15] 

W{n,T)=p c (vl-pT p f. (11) 

Inserting equation (11) into relation (8) and using equation (9), gives the 
relation between fa and 9 



fa = 2/3T p ^2p c K sign(fl - -)^cos- (l - cos - J (12) 

where a = cos _1 (sin 2 6) and the function sign returns the sign of its argument. 

We impose a no-slip boundary condition on the velocity. Because a collision 
takes place on the boundary the usual bounce-back condition must be extended 
to ensure mass conservation (see [11] for a wider discussion). This is done by a 
suitable choice of the rest field, /o, to correctly balance the mass of the system. 



3 Spreading on a homogeneous surface 



We consider a 80 x 80 x 40 lattice on which a spherical drop of radius Rq = 16 
just touches a flat surface at z = 0. Unless otherwise specified the temperature 
is T = 0.4 which leads to two phases of density m = 4.128 and n g — 2.913. 
Fig. 1 shows how the droplet evolves in time to reach an equilibrium shape with 
contact angle 60°. 
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Fig. 1. Spreading of a spherical droplet of radius Ro = 16 on a 80 x 80 x 40 lattice. 
The equilibrium contact angle is 60°. 

We first present a check on the accuracy of the equilibrium properties of 
the model. Fig. 2 reports a comparison between two methods of measuring the 
contact angle. 9 V is the contact angle obtained from equation (9) with the surface 
tensions measured at equilibrium. 8 9 is the contact angle measured from the 
profile of the simulated droplet once equilibrium is reached. The agreement is 
good. Small errors results from the difficulty of a direct measurement of the 
contact angle on a discrete lattice. 

The shape of the area formed by the contact of a droplet with a homogeneous 
substrate is a disk. Its radius R c is a quantity which is rather simple to measure 
and has consequently attracted the attention of many scientists, see [17] and 
references therein. The time evolution of R c has been found to follow a power 
law R c — mt n ' 2 . The exponent n has been widely reported in the literature but 
with no consistent result. Marmur [17] in his review lists exponents between 0.06 
and 0.6. The value of m appears to be related to the droplet volume. 

Fig. 3 shows the time evolution of R c for different values of the viscosity and 
the surface tension. The curves correspond to a value n — 0.56 which is within 
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Fig. 2. Comparison between equilibrium contact angles 8 y and 6 3 (denned in the text) 
on a 110 x 110 x 50 lattice. The input contact angles are set from 30° to 140° every 
10°. t = 1.0 and k = 0.003. The initial droplet has a radius Ro = 18. 80 000 iterations 
were used to reach each equilibrium. The dashed line is the expected agreement. 



the range reported in the literature. The power law is independent of the surface 
tension and the viscosity. 

Indeed if the evolution curves are plotted as a function of the dimensionless 
time [18] 

t^ t * = ^-t (13) 

the data collapses onto a single curve as shown in fig. 3(b). Experimental data 
taken from [18] shows similar behaviour. 



4 Spreading on heterogeneous surfaces 

Almost any surface will contain physical and chemical inhomogcncities which will 
affect the spreading of a mesoscopic droplet. It has recently become feasible to 
fabricate surfaces with well-defined chemical properties on micron length scales 
and it is becoming possible to perform well-controlled experiments which probe 
the behaviour of mesoscopic droplets on chemically and physically heterogeneous 
substrates. Thus it is particularly interesting to develop techniques to model the 
effect of these surfaces on the spreading properties of a droplet. 

One of the simplest heterogeneous surfaces can be formed by alternating 
stripes of materials with different wetting properties. The static properties of 
droplets on such substrates have been discussed [19-21]. However less attention 
has been paid to the dynamics of the spreading. 

In this section we consider heterogeneous surfaces formed by alternating hy- 
drophilic and hydrophobic stripes. They are characterised by widths w p hi,w p ho 
and contact angles 9 p hi, p h o respectively. Fig. 4 presents the behaviour of a 




Fig. 3. Time evolution of the radius of the droplet base R c on a 90 x 90 x 50 lattice (a) 
as a function of time t (b) as a function of dimensionless time t* . The contact angle is 
set to 60° and Ro = 16. The solid line is the result of laboratory experiments [18]. 



three-dimensional droplet spreading on such a surface with Q p hi = 50°, p h o = 
110°, Wphi — 6, Wpho = 5. The droplet has an initial radius Rq — 18. 

It is apparent from the figure that the behaviour of the droplet depends on 
whether it is on a hydrophobic or a hydrophilic stripe. The equilibrium shape of 
the contact line shown in fig. 4(b) reflects the pattern of the underlying substrate 
which is comparable to that found in laboratory experiments [20]. 

The time evolution of the contact line is also shown in fig. 4(b). Note that 
its velocity decreases smoothly in the y-direction parallel to the stripes but 
not in the x-direction where it moves faster on the hydrophilic than on the 
hydrophobic stripes. Note also that the droplet remains symmetric about an 
axis perpendicular to the stripes but that the shape becomes asymmetric about 
an axis parallel to the stripes, depending upon the initial position of the center 
of the droplet. 

Observation of the movement of the contact line in the x-direction shows that 
in a hydrophilic region the contact angle tends to decrease and the velocity of the 
contact line increase. When the contact line reaches the boundary its progress 
is stopped and the contact angle increases until it is large enough to cross the 
hydrophobic stripe. 

It has been proposed that an equilibrium droplet on such a surface has a 
spatially averaged contact angle following Cassie's law [21] 

cos e = p phl cos e phi + p pho cos e pho (14) 

where p p hi and p p ho are the proportion of the substrate area which are hy- 
drophilic or hydrophobic respectively. However, this relation is not universally 
accepted. In particular it has been argued that there should be a dependence on 
the relative size of the droplet and the surface stripes [19]. Fig. 4(c)-(d) show 



characteristic angles for the droplet considered here. Their average is 76.5° which 
is close to the one predicted by Cassie's law, 78.7°. 

5 Conclusion 

We have used a three-dimensional lattice Boltzmann algorithm to model the 
spreading of a mesoscopic droplet. By incorporating the Cahn theory of wetting 
into the simulation we obtain a way of easily tuning the contact angle of the 
droplet on the substrate. This gives us the ability to simulate spreading on both 
homogeneous and heterogeneous surfaces. 

The approach provides a well-controlled way of investigating the dependence 
of the spreading on such properties as the droplet volume, contact angles and 
the substrate geometry. Further work is in progress to systematically determine 
how these parameters affect the velocity and shapes of the spreading droplets. 
It would also be interesting to investigate the effect of physical inhomogeneities 
on the spreading and to consider a droplet spreading on a porous surface. 

A particular aim of the work will be to compare the results to forthcoming 
experiments on substrates which have chemical patterning on mesoscopic length 
scale. This will allow increased understanding of both the experimental results 
and the model assumptions. For example we assume, as is the standard practice, 
no-slip boundary conditions on the velocity. These may not be appropriate on 
short length scales near a contact line. Moreover the liquid-gas density differ- 
ence in lattice Boltzmann models is very small compared to real fluids and it is 
important to undertake further work to assess the effect of this on the modelling 
results. 
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